% This program can be used to compute the value of DeltaTilde, defined
% as the maximum magnitude of a Ponzi scheme at time 0 that is consistent
% with the economy being depressed while converging to the secular
% stagnation steady state.
clear;
format long;

T=32; %Number of years of horizon

%% Calibration

%Exogenous parameters
rho=0.04; %Discount rate
n=-0.007; %Population growth rate
g=0.013; %Rate of technical progress
frisch=0.85; %Frisch elasticity of labor supply
piW=g; %Downward wage rigidity

%Targets to match
LaborSupply=1; %Labor supply in (neoclassical) steady state
ConsumptionGap=0.1; % -(cSS-cFB)/cFB
MaxPonzi=1.25; %Ratio of Ponzi scheme to GDP in the Ponzi steady state
HalfLife=2; %Years necessary for the inflation anchor to close half the gap with the (constant) rate of inflation
PhillipsSlope=-0.3; %Elasticity of nominal wage inflation with respect to unemployment
MinWealth=2; %(Negative of) minimum wealth relative to consumption under secular stagnation

%Endogenous parameters
kL=LaborSupply^(-1-1/frisch); %Scale parameter of labor supply function
theta=log(2)/HalfLife; %Speed of inflation deanchoring
beta=-PhillipsSlope*frisch; %Slope coefficient of the Phillips curve

%Neoclassical steady state & Secular stagnation steady state
LNeo=LaborSupply;
LSS=(1-ConsumptionGap)*LNeo;
NaturalRate=g-rho*(LNeo/LSS-1)-LNeo/LSS*piW;
lSS=(kL*LSS)^(-frisch); %Labor supply

%Ponzi steady state
deltaPonzi=MaxPonzi*LNeo;

%Calibration of the preference for wealth
RefWealth=-MinWealth*LSS;
sigma=log(LNeo*(rho+piW)/(LSS*(rho-n)))/log(1-deltaPonzi/RefWealth);
kW=((rho-n)/LNeo)*(deltaPonzi-RefWealth)^sigma;

%% Solving for the value of delta tilde

DeTilde=0.560794; %Guess a value of delta tilde

func=@(t) kW*((exp(-(piW+n)*t)*DeTilde-RefWealth).^(-sigma)).*exp(-(piW+rho)*t);

integral(func,0,Inf)-kW*(deltaPonzi-RefWealth)^(-sigma)/(rho-n) %If the guess is correct, this should be equal to zero (since this is the equation that defines delta tilde)
